load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;**************************************
begin

season = getenv("season")
varname = getenv("varname")
iregion = toint(getenv("iregion"))

season = "DJF"
varname = "pr5d"
iregion = 5


;----------------------------
Ntheta = 100
if (varname .eq. "pr1d") then
  theta_scale_theory := fspan(0.2,20,Ntheta)
else if (varname .eq. "pr5d") then
  theta_scale_theory := fspan(0.1,10,Ntheta)
else if (varname .eq. "pr30d") then
  theta_scale_theory := fspan(0.03,3,Ntheta)
end if
end if
end if

Nmu = 100
;mu_theory := fspan(0.8,1.493,Nmu)
mu_theory := fspan(0.95,1.643,Nmu)



;---------------------

Nmodel = 60
N1 = 35        ;CMIP5
N2 = 25        ;CMIP6


Nmodel_valid = 52           ;reach 3C warming
Nmodel_cri = toint(Nmodel_valid*0.80)

nlat = 180
nlon = 360

Nlevel = 6
warming_levels = (/"0C","0.5degC","1degC","2degC","3degC","4degC"/)

level_idx := 4        ;3C
index_idx := 1        ;R95p


;*******************************************
;regional
;*******************************************
lonW := (/0,0,80,115,0,260/)
lonE := (/360,360,150,180,25,300/)
latS := (/45,-70,45,30,50,50/)
latN := (/70,-45,70,45,70,70/)

subregions := (/"NH_extratropics","SH_extratropics","Northern_Asia","WNP_monsoon","Europe","Eastern_NAmerica"/)
subregions_plot := (/"NH mid-latitudes","SH mid-latitudes","Northern Asia","West North Pacific","Europe","Eastern North America"/)
domain := (/"ocean","ocean","land","ocean","land","land"/)
Nregion := dimsizes(lonW)



;*************************************************
;--------------read gamma_parameter---------------
diri := "/WORK/zhangwx/EConstraint/data/"
k_shape := new((/Nmodel,Nlevel,nlat,nlon/),float)
theta_scale := k_shape

f1 := addfile(diri+"CMIP5/pr/"+varname+"/"+season+"/gamfit_wet0.1/multimodel/gamfit_parameter_"+varname+"_multimodel_warming_levels_"+season+".nc","r")
k_shape(0:N1-1,:,:,:) = f1->k_shape              ;(35model, 6level, lat, lon)
theta_scale(0:N1-1,:,:,:) = f1->theta_scale

f2 := addfile(diri+"CMIP6/pr/"+varname+"/"+season+"/gamfit_wet0.1/multimodel/gamfit_parameter_"+varname+"_multimodel_warming_levels_"+season+".nc","r")
k_shape(N1:N1+N2-1,:,:,:) = f2->k_shape              ;(25model, 6level, lat, lon)
theta_scale(N1:N1+N2-1,:,:,:) = f2->theta_scale


;---------------------------------------
mean := k_shape*theta_scale               ;(model level lat lon)
copy_VarCoords(k_shape,mean)

mean_control := mean(:,0,:,:)             ;(model lat lon)      mean_pd
mean_control_MME := dim_median_n(mean_control,0)         ;(lat lon)
copy_VarCoords(mean_control(0,:,:),mean_control_MME)

;---------------------------------

alpha := mean(:,level_idx,:,:)/mean(:,0,:,:)        ;Mean1/Mean0   ;(model lat lon)
copy_VarCoords(mean(:,0,:,:),alpha)

alpha_MME := dim_median_n(alpha,0)                  ;(lat lon)
copy_VarCoords(alpha(0,:,:),alpha_MME)


;***************************************
;mask to land
;***************************************
a := addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")
lsdata := a->LSMASK

lsm := landsea_mask(lsdata,k_shape&lat,k_shape&lon)
lsm1 := conform(k_shape,lsm,(/2,3/))                 ;(/Nmodel,Nlevel,lat,lon/)

mean_control_MME_land := where(((lsm1(0,0,:,:) .eq. 1) .or. (lsm1(0,0,:,:) .eq. 3)), mean_control_MME, mean_control_MME@_FillValue)
copy_VarCoords(mean_control_MME,mean_control_MME_land)

alpha_MME_land := where(((lsm1(0,0,:,:) .eq. 1) .or. (lsm1(0,0,:,:) .eq. 3)), alpha_MME, alpha_MME@_FillValue)
copy_VarCoords(alpha_MME,alpha_MME_land)


;***********************************************
;regional mean
;***********************************************
lat := k_shape&lat
rad = 4*atan(1.0)/180.
clat := cos(lat*rad)
clat!0 = "lat"
clat&lat = lat

if domain(iregion).eq."ocean" then
  mean_control_MME_regional := wgt_areaave_Wrap(mean_control_MME({latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  alpha_MME_regional        := wgt_areaave_Wrap(       alpha_MME({latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
else if domain(iregion).eq."land" then
  mean_control_MME_regional := wgt_areaave_Wrap(mean_control_MME_land({latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  alpha_MME_regional        := wgt_areaave_Wrap(       alpha_MME_land({latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
end if
end if



;***************************************************************
;cal PR in nu & scale space based on MME mean_control and alpha for this region
;***************************************************************
prob0 = 0.05             ;R95

PR := new((/Nmu,Ntheta/),float)
PR = PR@_FillValue

;------------------------
  do itheta = 0,Ntheta-1
    do imu = 0,Nmu-1

      mean_0 := mean_control_MME_regional
      theta_scale_0 := theta_scale_theory(itheta)
      k_shape_0 := mean_0/theta_scale_0                ;mean=k*theta

      mu_tmp := mu_theory(imu)
      delta := (alpha_MME_regional/mu_tmp-1)*k_shape_0
      k_shape_1 := k_shape_0+delta
      theta_scale_1 := theta_scale_0*mu_tmp

      if (k_shape_1 .le. 0) then                 ;both shape & scale must >0
        continue
      end if

      RX_threshold_0 := cdfgam_x(1-prob0, k_shape_0, 1./theta_scale_0)           ;event threshold at current climate
      prob1 := 1.0 - cdfgam_p(RX_threshold_0, k_shape_1, 1./theta_scale_1)

      PR(imu,itheta) = prob1/prob0

    end do
  end do

;-------------------------
PR!0 = "mu"
PR&mu = mu_theory
PR!1 = "scale"
PR&scale = theta_scale_theory


;***************************************************************
;model data in parameter space -- dots
;***************************************************************

theta_scale(:,0,:,:) = where( theta_scale(:,0,:,:) .eq. 0, theta_scale@_FillValue, theta_scale(:,0,:,:) )

mu_model := theta_scale(:,level_idx,:,:)/theta_scale(:,0,:,:)               ;theta1/theta0   (model lat lon)
copy_VarCoords(theta_scale(:,0,:,:),mu_model)


;******************mask to land********************

mu_model_land := where(((lsm1(:,0,:,:) .eq. 1) .or. (lsm1(:,0,:,:) .eq. 3)), mu_model, mu_model@_FillValue)
copy_VarCoords(mu_model,mu_model_land)

theta_scale_land := where(((lsm1(:,:,:,:) .eq. 1) .or. (lsm1(:,:,:,:) .eq. 3)), theta_scale, theta_scale@_FillValue)
copy_VarCoords(theta_scale,theta_scale_land)

;--------------------------------------------
if domain(iregion).eq."ocean" then
  mu_model_regional          := wgt_areaave_Wrap(     mu_model(:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  theta_scale_model_regional := wgt_areaave_Wrap(theta_scale(:,0,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
else if domain(iregion).eq."land" then
  mu_model_regional          := wgt_areaave_Wrap(     mu_model_land(:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
  theta_scale_model_regional := wgt_areaave_Wrap(theta_scale_land(:,0,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)
end if
end if


;***************************************************************
;PR estimates from direct model output
;***************************************************************

change_frq_model := new((/Nmodel,nlat,nlon/),float)

f1 := addfile(diri+"CMIP5/pr/"+varname+"/"+season+"/RX_pct/multimodel/change_frequency_RX_pct_"+varname+"_multimodel_warming_levels_"+season+".nc","r")
change_frq_model(0:N1-1,:,:) = f1->change_frq_RX_pct(:,level_idx,index_idx,:,:)               ;(35model, lat, lon)

f2 := addfile(diri+"CMIP6/pr/"+varname+"/"+season+"/RX_pct/multimodel/change_frequency_RX_pct_"+varname+"_multimodel_warming_levels_"+season+".nc","r")
change_frq_model(N1:N1+N2-1,:,:) = f2->change_frq_RX_pct(:,level_idx,index_idx,:,:)        ;(25model lat lon)

PR_model := change_frq_model/100.+1.
copy_VarCoords(change_frq_model,PR_model)

;******************mask to land********************

PR_model_land := where(((lsm1(:,0,:,:) .eq. 1) .or. (lsm1(:,0,:,:) .eq. 3)), PR_model, PR_model@_FillValue)
copy_VarCoords(PR_model,PR_model_land)

;--------------------------------------
if domain(iregion).eq."ocean" then
  PR_model_regional := wgt_areaave_Wrap(PR_model(:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)       ;{model}
else if domain(iregion).eq."land" then
  PR_model_regional := wgt_areaave_Wrap(PR_model_land(:,{latS(iregion):latN(iregion)},{lonW(iregion):lonE(iregion)}),clat({latS(iregion):latN(iregion)}),1.0,0)       ;{model}
end if
end if



;*****************************************************
wks := gsn_open_wks("pdf","/WORK/zhangwx/EConstraint/pic/gamma_test/ModelFit/Theory_Model/onlyCMIP/round1/PR_R95_gamma_distribution_transform2_"+varname+"_"+season+"_"+subregions(iregion))

;gsn_define_colormap(wks,"WhiteBlueGreenYellowRed")
;cmap := read_colormap_file("WhiteBlueGreenYellowRed")
cmap := read_colormap_file("MPL_YlGnBu")

plot_pm := new(Nmodel,graphic)
plot_pm_hollow := new(Nmodel,graphic)


res = True
res@gsnFrame = False
res@gsnDraw = False

res@cnFillOn = True
res@cnLinesOn = False
res@cnLineLabelsOn = False
res@cnInfoLabelOn = False

;res@mpProjection = "Robinson"
;res@mpPerimOn = False

;res@gsnSpreadColors = True
;res@gsnSpreadColorStart = 63
;res@gsnSpreadColorEnd = 4
res@cnMissingValFillColor = "grey"

res@tmXTOn = False
res@tmYROn = False

  res@tmBorderThicknessF = 1.0
  res@tmXBMajorThicknessF = 1.0
  res@tmXBMinorThicknessF = 1.0
  res@tmYLMajorThicknessF = 1.0
  res@tmYLMinorThicknessF = 1.0

  res@tmXBLabelFontHeightF = .017
  res@tmYLLabelFontHeightF = .017

res@gsnLeftStringFontHeightF = 0.021
res@gsnRightStringFontHeightF = 0.021


;res@cnFillPalette = "BlueWhiteOrangeRed"
res@lbLabelBarOn = False
;res@pmLabelBarOrthogonalPosF = -0.03
;res@pmLabelBarHeightF = 0.06

;res@tiMainFontHeightF = 0.023
;res@tiMainFontThicknessF = 2

res@gsnLeftString = ""
res@gsnRightString = ""

res@tiYAxisString = "~F33~n"       ;nu
res@tiXAxisString = "~F33~q"       ;theta_scale
res@tmYLPrecision = 3

;res@trYMinF = 0.6
if (max(theta_scale_model_regional)*1.5) .lt. max(theta_scale_theory) then
  res@trXMaxF = (max(theta_scale_model_regional)*1.5)
end if

N_PR_bin = 41
;if (varname .eq. "pr1d") then
;  PR_bin := fspan(1,1.6,N_PR_bin)
;else if (varname .eq. "pr5d") then
;  PR_bin := fspan(1,2.,N_PR_bin)
;else if (varname .eq. "pr30d") then
;  PR_bin := fspan(1,3.,N_PR_bin)
;end if
;end if
;end if
if (subregions(iregion).eq."NH_extratropics") then
  PR_bin := fspan(1,1.7,N_PR_bin)
else if (subregions(iregion).eq."SH_extratropics") then
  PR_bin := fspan(1.3,2.1,N_PR_bin)
else if (subregions(iregion).eq."Northern_Asia") then
  PR_bin := fspan(1,1.8,N_PR_bin)
else if (subregions(iregion).eq."WNP_monsoon") then
  PR_bin := fspan(1.2,2.,N_PR_bin)
else if (subregions(iregion).eq."Europe") then
  PR_bin := fspan(1.2,2.5,N_PR_bin)
else if (subregions(iregion).eq."Eastern_NAmerica") then
  PR_bin := fspan(1.64,2.84,N_PR_bin)
end if
end if
end if
end if
end if
end if


;--------------for southern_oean--------------
;if (varname .eq. "pr1d").and.(iregion .eq. 2)then      ;;; for southern_oean
;   PR_bin := fspan(1.16,1.56,N_PR_bin)
;end if
;if (varname .eq. "pr5d").and.(iregion .eq. 2)then      ;;; for southern_oean
;   PR_bin := fspan(1.36,2.,N_PR_bin)
;end if
;if (varname .eq. "pr30d").and.(iregion .eq. 2)then      ;;; for southern_oean
;   PR_bin := fspan(1.5,4,N_PR_bin)
;end if


;color_bin := ispan(2,248,6)            ;42 colors
color_bin := ispan(2,127,3)            ;42 colors

res@cnLevelSelectionMode = "ExplicitLevels"
res@cnLevels := PR_bin
res@cnFillColors = cmap(color_bin,:)

;-------------------------------------------------
leftst = (/"a","b","c","d","e","f"/)

res@gsnLeftString = subregions_plot(iregion)
res@gsnRightString = "Mean="+sprintf("%4.2f",mean_control_MME_regional)+", ~F33~a~F21~="+sprintf("%4.2f",alpha_MME_regional)
plot = gsn_csm_contour(wks,PR,res)


txres = True
txres@txFont = "helvetica-bold"
txres@txFontHeightF = 0.021
gsn_text_ndc(wks,leftst(iregion),0.084,0.86,txres)


;*******************************************************
;model data in parameter space
;*******************************************************


pmres = True
pmres@gsMarkerIndex =  16
pmres@gsMarkerSizeF = 0.01

pmres_hollow = True
pmres_hollow@gsMarkerIndex =  4
pmres_hollow@gsMarkerSizeF = 0.01
pmres_hollow@gsMarkerThicknessF = 0.5


do imodel = 0,Nmodel-1
  if (.not.(ismissing(PR_model_regional(imodel)))) then

    bin_idx := min( ind( PR_model_regional(imodel) .le. PR_bin ) )
    if ( PR_model_regional(imodel) .ge. PR_bin(N_PR_bin-1) ) then
      bin_idx := N_PR_bin-1
    end if
    ;print(PR_bin(bin_idx)+"  "+PR_model_regional(imodel)+"  "+bin_idx)
    pmres@gsMarkerColor = cmap(color_bin(bin_idx),:)
    plot_pm(imodel) = gsn_add_polymarker(wks,plot,theta_scale_model_regional(imodel),mu_model_regional(imodel),pmres)

    plot_pm_hollow(imodel) = gsn_add_polymarker(wks,plot,theta_scale_model_regional(imodel),mu_model_regional(imodel),pmres_hollow)

  end if
end do





;----------------------------------
;--------------------
resP = True
;resP@gsnPanelTop = 0.95
;resP@gsnPanelLeft = 0.05
resP@gsnPanelLabelBar  = True
resP@gsnPanelYWhiteSpacePercent = 3.
resP@gsnPanelXWhiteSpacePercent = 4.
resP@lbOrientation = "Vertical"
resP@lbLabelFontHeightF = 0.018
resP@lbLabelStride = 4
resP@pmLabelBarHeightF = 0.6
resP@pmLabelBarWidthF = 0.1
resP@pmLabelBarOrthogonalPosF = 0.02
resP@pmLabelBarParallelPosF = 0.0

if iregion .eq. 3 then
  resP@pmLabelBarHeightF = 0.59
  resP@pmLabelBarOrthogonalPosF = -0.004
  resP@pmLabelBarParallelPosF = 0.0
end if

gsn_panel(wks,plot,(/1,1/),resP)


end
